For this analysis, we used the following parameters:

pars <- as.data.frame(do.call("rbind", params))
colnames(pars) <- "parameters_used"
pars

Load Seurat object.

git_path <- params$git_path
figure_out_path <- params$figure_out_path

sobj <- readRDS(file=file.path(git_path,params$object))

Load DESeq2 pseudobulk results. For full documentation of DESeq2 analysis, refer to Scritps/DESeq2.Rmd.

res_df <- read.delim(file.path(git_path,"Figures/data/de_results_pseudoid_analysis.tsv"))


selected_contrasts <- c("12w_nonSPF_5days_fx_vs_12w_nonSPF_2days_fx",
                        "26w_nonSPF_5days_fx_vs_26w_nonSPF_2days_fx",
                        "26w_SPF_5days_fx_vs_26w_SPF_2days_fx",
                        "26w_nonSPF_2days_fx_vs_12w_nonSPF_2days_fx",
                        "26w_SPF_2days_fx_vs_12w_nonSPF_2days_fx",
                        "26w_SPF_2days_fx_vs_26w_nonSPF_2days_fx",
                        "12w_nonSPF_5days_distal_vs_12w_nonSPF_5days_fx",
                        "12w_nonSPF_5days_proximal_vs_12w_nonSPF_5days_fx",
                        "12w_nonSPF_5days_distal_vs_12w_nonSPF_5days_proximal",
                        "26w_nonSPF_5days_distal_vs_26w_nonSPF_5days_fx",
                        "26w_nonSPF_5days_proximal_vs_26w_nonSPF_5days_fx",
                        "26w_nonSPF_5days_distal_vs_26w_nonSPF_5days_proximal",
                        "26w_SPF_5days_distal_vs_26w_SPF_5days_fx",
                        "26w_SPF_5days_proximal_vs_26w_SPF_5days_fx",
                        "26w_SPF_5days_distal_vs_26w_SPF_5days_proximal")

cs <- data.frame(contrast=unique(res_df$contrast),
                 new=selected_contrasts)

res_df <- res_df %>%
  left_join(cs) %>%
  dplyr::rename(old_contrast=contrast,
         contrast=new)

gene_info <- read.delim(file.path(git_path,params$gene_info), header=F)
colnames(gene_info) <- c("gene_id","gene_name","gene_type")

res_df <- res_df %>%
  left_join(gene_info)
for (comparison in unique(res_df$contrast)){
  p <- res_df %>%                  
    dplyr::filter(!is.na(padj), contrast == comparison) %>%
    ggplot(aes(x=log2FoldChange,y=-log10(pvalue),color=padj < .05)) + 
    geom_point(size=.5, shape=1) +
    geom_label_repel(data=res_df %>%
                       dplyr::filter(contrast == comparison,padj < .05),
                     aes(label=gene_name),
                     segment.size=.25,segment.alpha=.5,size=3,color='black',max.overlaps=40,
                     label.padding=0.05,label.size=NA) +
    scale_color_manual(values=c('black','red')) +
    theme_classic() +
    theme(legend.position='none',
          strip.background=element_blank())+
    ggtitle(comparison)+
    facet_wrap(~cell_type, ncol = 2, scales = "free")
  print(p)
  
}

res_df %>%
  dplyr::arrange(padj) %>%
  dplyr::filter(!is.na(padj) & (padj < .05)) %>%
  dplyr::mutate_if(is.numeric, signif, 2) %>%
  DT::datatable(rownames=FALSE, extensions = 'Buttons', options = list(dom='frtBip', buttons = c('csv')))
library(tmod)

# this is to document how the mouse orthologues were mapped for the tmod pathways

# data(tmod)
# 
# library(orthomapper)
# mousehum <- orthologs(10090, 9606)
# mousehum$h_symbol <- entrez_annotate(mousehum$T.9606)$SYMBOL
# mousehum$m_symbol <- entrez_annotate(mousehum$T.10090)$SYMBOL
# 
# mousehum <- mousehum[order(mousehum$h_symbol),]
# rows_tmod <- which(mousehum$h_symbol %in% tmod$gv)
# mgenes <- mousehum$m_symbol[rows_tmod]
# 
# list_tmod <- sapply(tmod$gs2gv, function(x) tmod$gv[x] )
# list_tmod <- sapply(list_tmod, function(x) x[order(x)])
# list_tmod <- sapply(list_tmod, function(x) x[which(x %in%  mousehum$h_symbol)])
# list_tmod <- sapply(list_tmod, function(x) mousehum[which(mousehum$h_symbol %in% x), "m_symbol"])
# list_tmod <- sapply(list_tmod, function(x) which(mgenes %in% x))
# 
# tmod_mm <- tmod
# 
# tmod_mm$gv <- mgenes
# tmod_mm$gs2gv <- list_tmod

# save tmod mouse object
# saveRDS(tmod_mm, "tmod_mm.rds")


tmod_mm <- readRDS(file.path(git_path,"Figures/data/tmod_mm.rds"))

df2tmod <- function(df, gene_id_col=ncol(df), module_id_col=1, module_title_col=2) {

  require(tmod, quietly=TRUE)

  gene_ids <- df[[gene_id_col]]
  module_ids <- df[[module_id_col]]
  m2g <- tapply(gene_ids, module_ids, list)
  df[[gene_id_col]] <- NULL
  df <- df[!duplicated(df[[module_id_col]]), ]
  colnames(df)[module_id_col] <- "ID"
  colnames(df)[module_title_col] <- "Title"

  makeTmod(modules=df, modules2genes=m2g)
}

df <- as.data.frame(msigdbr::msigdbr(species='Mus musculus'))
df <- df[ , c("gs_name", "gs_id", "gs_cat", "gs_subcat", "gene_symbol") ]
colnames(df) <- c("Title", "ID", "Category", "Subcategory", "GeneID")
msig <- df2tmod(df[!is.na(df$GeneID),], gene_id_col=ncol(df), module_id_col=2, module_title_col=1)
sel <- (msig$gs$Category %in% c("H"))  | (msig$gs$Subcategory %in% c('CP:REACTOME','GO:BP','CP:KEGG'))
res <- res_df %>%
  filter(gene_type=="protein_coding")
res <- split(res, f = res$contrast)

genes <- lapply(res, function(x) x %>%
                  filter(!is.na(padj)) %>%
                  arrange(pvalue))
genes_reduced <- genes["26w_nonSPF_5days_fx_vs_26w_nonSPF_2days_fx"]

lgenes <- lapply(genes_reduced, function(x) x$gene_name)

res.tmod <- lapply(lgenes, tmodCERNOtest, mset = tmod_mm)

res.tmod.filtered <- lapply(res.tmod, function(x) subset(x, adj.P.Val < 1e-3 & AUC > .75))

Figure 3A

#################### Volcano plot #######################
DF <- res_df[grepl("26w_nonSPF_5days_fx_vs_26w_nonSPF_2days_fx", res_df$contrast), ]
# add a column of NAs
DF$direction <- "N"
# if log2Foldchange > 0.5 and pvalue < 0.05, set as "UP" 
DF$direction[(DF$log2FoldChange > 0.5 & DF$padj < 0.1)] <- "up"
DF$direction[(DF$log2FoldChange  < (-0.5) & DF$padj < 0.1)] <- "down"

(p12 <- DF %>% 
    dplyr::filter(gene_type=="protein_coding") %>%
    dplyr::filter(!is.na(padj)) %>%
    ggplot(aes(x=log2FoldChange,y=-log10(pvalue),col=direction)) + 
    geom_point(size=.5, shape=1) +
    geom_label_repel(data=DF %>%
                       dplyr::filter(gene_type=="protein_coding") %>% 
                       dplyr::filter(direction=="up"|direction=="down") ,
                     aes(label=gene_name),
                     segment.size=.25,segment.alpha=.5,size=2.5,color='black',max.overlaps=30,
                     label.padding=0.05,label.size=NA) +
    scale_color_manual(values=c("blue","grey","red"))+ theme_classic() + 
    theme(text = element_text(family = "Arial"),
        plot.title = element_text(size = 8,hjust = 0.5,),
                     axis.title = element_text(size = 8),  
                     axis.text = element_text(size = 8),
                     axis.ticks.x=element_blank(),
                     legend.text = element_text(size =8),
                     legend.title = element_text(size = 8),
                     strip.text.x = element_blank()) + 
    labs(title = "",x = "log2FoldChange (5 days vs. 2 days)") +
  ggtitle("Immune-aged: hematoma"))

Figure S3A

#################### Volcano plot #######################
DF <- res_df[grepl("26w_SPF_5days_fx_vs_26w_SPF_2days_fx", res_df$contrast), ]
# add a column of NAs
DF$direction <- "N"
# if log2Foldchange > 0.5 and pvalue < 0.05, set as "UP" 
DF$direction[(DF$log2FoldChange > 0.5 & DF$padj < 0.1)] <- "up"
DF$direction[(DF$log2FoldChange  < (-0.5) & DF$padj < 0.1)] <- "down"

(sf3a <- DF %>% 
    dplyr::filter(gene_type=="protein_coding") %>%
    dplyr::filter(!is.na(padj)) %>%
    ggplot(aes(x=log2FoldChange,y=-log10(pvalue),col=direction)) + 
    geom_point(size=.5, shape=1) +
    geom_label_repel(data=DF %>%
                       dplyr::filter(gene_type=="protein_coding") %>% 
                       dplyr::filter(direction=="up"|direction=="down") ,
                     aes(label=gene_name),
                     segment.size=.25,segment.alpha=.5,size=2.5,color='black',max.overlaps=30,
                     label.padding=0.05,label.size=NA) +
    scale_color_manual(values=c("blue","grey","red"))+ theme_classic() + 
    theme(text = element_text(family = "Arial"),
        plot.title = element_text(size = 8,hjust = 0.5,),
                     axis.title = element_text(size = 8),  
                     axis.text = element_text(size = 8),
                     axis.ticks.x=element_blank(),
                     legend.text = element_text(size =8),
                     legend.title = element_text(size = 8),
                     strip.text.x = element_blank()) + 
    labs(title = "",x = "log2FoldChange (5 days vs. 2 days)")+
  ggtitle("Aged: hematoma"))

Figure S3B

#################### Volcano plot #######################
DF <- res_df[grepl("12w_nonSPF_5days_fx_vs_12w_nonSPF_2days_fx", res_df$contrast), ]
# add a column of NAs
DF$direction <- "N"
# if log2Foldchange > 0.5 and pvalue < 0.05, set as "UP" 
DF$direction[(DF$log2FoldChange > 0.5 & DF$padj < 0.1)] <- "up"
DF$direction[(DF$log2FoldChange  < (-0.5) & DF$padj < 0.1)] <- "down"

(sf3b <- DF %>% 
    dplyr::filter(gene_type=="protein_coding") %>%
    dplyr::filter(!is.na(padj)) %>%
    ggplot(aes(x=log2FoldChange,y=-log10(pvalue),col=direction)) + 
    geom_point(size=.5, shape=1) +
    geom_label_repel(data=DF %>%
                       dplyr::filter(gene_type=="protein_coding") %>% 
                       dplyr::filter(direction=="up"|direction=="down") ,
                     aes(label=gene_name),
                     segment.size=.25,segment.alpha=.5,size=2.5,color='black',max.overlaps=30,
                     label.padding=0.05,label.size=NA) +
    scale_color_manual(values=c("blue","grey","red"))+ theme_classic() + 
    theme(text = element_text(family = "Arial"),
        plot.title = element_text(size = 8,hjust = 0.5,),
                     axis.title = element_text(size = 8),  
                     axis.text = element_text(size = 8),
                     axis.ticks.x=element_blank(),
                     legend.text = element_text(size =8),
                     legend.title = element_text(size = 8),
                     strip.text.x = element_blank()) + 
    labs(title = "",x = "log2FoldChange (5 days vs. 2 days)") +
  ggtitle("Young: hematoma"))

Panel plot for the immune-aged hematoma.

res <- res_df %>%
  filter(gene_type=="protein_coding")
res <- split(res, f = res$contrast)

genes <- lapply(res, function(x) x %>%
                  filter(!is.na(padj)) %>%
                  arrange(pvalue))
genes_reduced <- genes["26w_nonSPF_5days_fx_vs_26w_nonSPF_2days_fx"]

lgenes <- lapply(genes_reduced, function(x) x$gene_name)

res.tmod <- lapply(lgenes, tmodCERNOtest, mset = tmod_mm)

res.tmod.filtered <- lapply(res.tmod, function(x) subset(x, adj.P.Val < 1e-3 & AUC > .75))


x <- genes_reduced$"26w_nonSPF_5days_fx_vs_26w_nonSPF_2days_fx"
sgenes <- tmodDecideTests(g = x$gene_name, lfc = x$log2FoldChange, pval = x$padj, lfc.thr = 0.5, pval.thr = 0.1, mset = tmod_mm)

names(sgenes) <- "26w_nonSPF_5days_fx_vs_26w_nonSPF_2days_fx"
p10 <- ggPanelplot(res.tmod.filtered, sgenes = sgenes,q_cutoff = 1e-7,cluster=F) + 
  theme(text = element_text(family = "Arial"),
        plot.title = element_text(size = 8,hjust = 0.5,),
        axis.title = element_text(size = 8),
        axis.text = element_text(size = 8),
        axis.ticks.x=element_blank(),
        legend.text = element_text(size =8),
        legend.title = element_text(size = 8),
        strip.text.x = element_blank()) + 
  labs(title = "Immune-aged: \nhematoma, \n5 days vs 2 days")
p10

For young and aged no enrichment is found, so no figures are generated here.

res <- res_df %>%
  filter(gene_type=="protein_coding")
res <- split(res, f = res$contrast)

genes <- lapply(res, function(x) x %>%
                  filter(!is.na(padj)) %>%
                  arrange(pvalue))
genes_reduced <- genes["12w_nonSPF_5days_distal_vs_12w_nonSPF_5days_fx"]

lgenes <- lapply(genes_reduced, function(x) x$gene_name)

res.tmod <- lapply(lgenes, tmodCERNOtest, mset = tmod_mm)

res.tmod.filtered <- lapply(res.tmod, function(x) subset(x, adj.P.Val < 1e-3 & AUC > .75))


x <- genes_reduced$"12w_nonSPF_5days_distal_vs_12w_nonSPF_5days_fx"
sgenes <- tmodDecideTests(g = x$gene_name, lfc = x$log2FoldChange, pval = x$padj, lfc.thr = 0.5, pval.thr = 0.1, mset = tmod_mm)

names(sgenes) <- "12w_nonSPF_5days_distal_vs_12w_nonSPF_5days_fx"
# sf3c <- ggPanelplot(res.tmod.filtered, sgenes = sgenes, q_cutoff = 1e-7,cluster=F) + 
#   theme(text = element_text(family = "Arial"),
#         plot.title = element_text(size = 8,hjust = 0.5,),
#         axis.title = element_text(size = 8),
#         axis.text = element_text(size = 8),
#         axis.ticks.x=element_blank(),
#         legend.text = element_text(size =8),
#         legend.title = element_text(size = 8),
#         strip.text.x = element_blank()) + 
#   labs(title = "Young: \nhematoma, \n5 days vs 2 days")
# sf3c
res <- res_df %>%
  filter(gene_type=="protein_coding")
res <- split(res, f = res$contrast)

genes <- lapply(res, function(x) x %>%
                  filter(!is.na(padj)) %>%
                  arrange(pvalue))
genes_reduced <- genes["26w_SPF_5days_fx_vs_26w_SPF_2days_fx"]

lgenes <- lapply(genes_reduced, function(x) x$gene_name)

res.tmod <- lapply(lgenes, tmodCERNOtest, mset = tmod_mm)

res.tmod.filtered <- lapply(res.tmod, function(x) subset(x, adj.P.Val < 1e-3 & AUC > .75))


x <- genes_reduced$"26w_SPF_5days_fx_vs_26w_SPF_2days_fx"
sgenes <- tmodDecideTests(g = x$gene_name, lfc = x$log2FoldChange, pval = x$padj, lfc.thr = 0.5, pval.thr = 0.1, mset = tmod_mm)

names(sgenes) <- "12w_nonSPF_5days_distal_vs_12w_nonSPF_5days_fx"
# sf3d <- ggPanelplot(res.tmod.filtered, sgenes = sgenes, q_cutoff = 1e-7,cluster=F) + 
#   theme(text = element_text(family = "Arial"),
#         plot.title = element_text(size = 8,hjust = 0.5,),
#         axis.title = element_text(size = 8),
#         axis.text = element_text(size = 8),
#         axis.ticks.x=element_blank(),
#         legend.text = element_text(size =8),
#         legend.title = element_text(size = 8),
#         strip.text.x = element_blank()) + 
#   labs(title = "Aged: \nhematoma, \n5 days vs 2 days")
# sf3d

Generating figure 4C

res <- res_df %>%
  filter(gene_type=="protein_coding")
res <- split(res, f = res$contrast)

genes <- lapply(res, function(x) x %>%
                  filter(!is.na(padj)) %>%
                  arrange(pvalue))
genes_reduced <- genes["26w_nonSPF_5days_fx_vs_26w_nonSPF_2days_fx"]

lgenes <- lapply(genes_reduced, function(x) x$gene_name)

res.tmod <- lapply(lgenes, tmodCERNOtest, mset = tmod_mm)

res.tmod.filtered <- lapply(res.tmod, function(x) subset(x, adj.P.Val < 1e-3 & AUC > .75))

v <- unlist(lgenes, use.names=FALSE)#2 for non SPF
down<-x %>% filter(log2FoldChange < 0.5, padj<0.1) %>% pull(gene_name)
modmetabo<-tmod_mm
modmetabo$gs[modmetabo$gs$ID %in% c("LI.M47.0"), ]
path_positions<- modmetabo$gs2gv[[which(modmetabo$gs$ID=="LI.M47.0")]]
genes<-modmetabo$gv[path_positions]

Genes_downs = down[down %in% genes]
Genes_all = v[v %in% genes]
gene.colors <- c("grey","grey","blue" ,"grey","blue" ,"grey", "blue",   "grey","grey","grey","grey", "grey","grey", "grey","grey","grey","grey", "grey", "grey", "grey","grey","grey","grey","grey","grey","grey" ,"grey","grey","grey","blue","grey" ,"grey","grey","grey","grey")
names(gene.colors)<- genes
######### Evidence plot ###############################
update_geom_defaults("text_repel", list(size = 2.5))
(p11 <- ggEvidencePlot(v, "LI.M47.0", mset = tmod_mm, gene.colors = gene.colors,filter=F) + 
    theme_classic() +  
    theme(text = element_text(family = "Arial",size=8),
          plot.title = element_text(size = 8,hjust = 0.5,),
          axis.title = element_text(size = 8), 
          axis.text = element_text(size = 8),
          axis.ticks.x=element_blank(),
          legend.text = element_text(size =8),
          legend.title = element_text(size = 8),
          legend.position = "none",
          strip.text.x = element_blank()) + labs(title = "Genes in LI.M47.0"))

Generating Figure E & F

############# Downsample ################################
# Setting seeds 
set.seed(10) 
Idents(sobj) <- "group"
Downsample<- subset(sobj, downsample=500)
############# Umap plot ################################

modmetabo<-tmod_mm
modmetabo$gs[modmetabo$gs$ID %in% c("LI.M47.0"), ]
path_positions<- modmetabo$gs2gv[[which(modmetabo$gs$ID=="LI.M47.0")]]
genes<-modmetabo$gv[path_positions]

Genes_Bcells= v[v %in% genes]

DefaultAssay(Downsample) <- "RNA"
Downsample <- AddModuleScore(Downsample,
                  features = list(Genes_Bcells),
                  name="Genes_Bcells")

Idents(Downsample)<-Downsample$timepoint
testdata <-subset(Downsample, idents = c("2days"))
Idents(testdata)<-testdata$compartments
test <-subset(testdata, idents = c("Hematoma"))
Idents(test)<-test$Age_raw
Immunoaged <- subset(test, idents = c("Immunoaged"))
Idents(Immunoaged)<-Immunoaged$CellTypes
# Plot module scores bcells 2 days 
p13<-FeaturePlot(Immunoaged,
                 features = "Genes_Bcells1",label = FALSE, repel = TRUE) +
  scale_colour_gradientn(colours = rev(brewer.pal(n = 11, name = "RdBu")), limits = c(0.0, 1.5)) +
  guides(colour=guide_colorbar(title="Gene module \nscore"))+
  theme(text = element_text(family = "Arial"),
        plot.title = element_text(size = 8),
        axis.title = element_text(size = 8),  
        axis.text = element_text(size = 8),
        legend.text = element_text(size =8),
        legend.title = element_text(size = 8),
        strip.text.x = element_text(size = 8),
        legend.position = "top",
        legend.justification = "center") + labs(title = "Immune-aged: hematoma, 2 days") & NoAxes()

p13

################################# 5 days ############################################
Idents(Downsample)<-Downsample$timepoint
testdata <-subset(Downsample, idents = c("5days"))
Idents(testdata)<-testdata$compartments
test <-subset(testdata, idents = c("Hematoma"))
Idents(test)<-test$Age_raw
Immunoaged <- subset(test, idents = c("Immunoaged"))
Idents(Immunoaged)<-Immunoaged$CellTypes
# Plot module scores bcells 5 days 
p14<-FeaturePlot(Immunoaged,
                 features = "Genes_Bcells1",label = FALSE, repel = TRUE) +
  scale_colour_gradientn(colours = rev(brewer.pal(n = 11, name = "RdBu")), limits = c(0.0, 1.5)) +
  guides(colour=guide_colorbar(title="Gene module \nscore"))+
  theme(text = element_text(family = "Arial"),
        plot.title = element_text(size = 8),
        axis.title = element_text(size = 8),  
        axis.text = element_text(size = 8),
        legend.text = element_text(size =8),
        legend.title = element_text(size = 8),
        strip.text.x = element_text(size = 8),
        legend.position = "top",
        legend.justification = "center") + labs(title = "Immune-aged: hematoma, 5 days")  & NoAxes()

p14

Interferon - which cell type? UMAP

# ############# Downsample ################################
# # Setting seeds 
# set.seed(10) 
# Idents(sobj) <- "group"
# Downsample<- subset(sobj, downsample=500)
# ############# Umap plot ################################
# 
# ifn <- read.delim(file.path(git_path,"Figures/data/signatures_mm_signatuR.txt"), header=T) %>%
#   filter(signature=="IfnResp") %>% 
#   pull(gene)
# 
# genes <- ifn
# 
# Genes_Bcells= v[v %in% genes]
# 
# DefaultAssay(Downsample) <- "RNA"
# Downsample <- AddModuleScore(Downsample,
#                   features = list(Genes_Bcells),
#                   name="IfnResponse")
# 
# Idents(Downsample)<-Downsample$timepoint
# testdata <-subset(Downsample, idents = c("2days"))
# Idents(testdata)<-testdata$compartments
# test <-subset(testdata, idents = c("Hematoma"))
# Idents(test)<-test$Age_raw
# Immunoaged <- subset(test, idents = c("Immunoaged"))
# Idents(Immunoaged)<-Immunoaged$CellTypes
# # Plot module scores bcells 2 days 
# p13<-FeaturePlot(Immunoaged,
#                  features = "IfnResponse1",label = FALSE, repel = TRUE) +
#   scale_colour_gradientn(colours = rev(brewer.pal(n = 11, name = "RdBu")), limits = c(0.0, 1.5)) +
#   guides(colour=guide_colorbar(title="Gene module \nscore"))+
#   theme(text = element_text(family = "Arial"),
#         plot.title = element_text(size = 8),
#         axis.title = element_text(size = 8),  
#         axis.text = element_text(size = 8),
#         legend.text = element_text(size =8),
#         legend.title = element_text(size = 8),
#         strip.text.x = element_text(size = 8),
#         legend.position = "top",
#         legend.justification = "center") + labs(title = "Immune-aged: hematoma, 2 days") & NoAxes()
# 
# p13
# 
# ################################# 5 days ############################################
# Idents(Downsample)<-Downsample$timepoint
# testdata <-subset(Downsample, idents = c("5days"))
# Idents(testdata)<-testdata$compartments
# test <-subset(testdata, idents = c("Hematoma"))
# Idents(test)<-test$Age_raw
# Immunoaged <- subset(test, idents = c("Immunoaged"))
# Idents(Immunoaged)<-Immunoaged$CellTypes
# # Plot module scores bcells 5 days 
# p14<-FeaturePlot(Immunoaged,
#                  features = "IfnResponse1",label = FALSE, repel = TRUE) +
#   scale_colour_gradientn(colours = rev(brewer.pal(n = 11, name = "RdBu")), limits = c(0.0, 1.5)) +
#   guides(colour=guide_colorbar(title="Gene module \nscore"))+
#   theme(text = element_text(family = "Arial"),
#         plot.title = element_text(size = 8),
#         axis.title = element_text(size = 8),  
#         axis.text = element_text(size = 8),
#         legend.text = element_text(size =8),
#         legend.title = element_text(size = 8),
#         strip.text.x = element_text(size = 8),
#         legend.position = "top",
#         legend.justification = "center") + labs(title = "Immune-aged: hematoma, 5 days")  & NoAxes()
# 
# p14

Generating Figure 4 D

library(ggpubr)
my_comparisons <- list(c("2days", "5days"))
DefaultAssay(Downsample) <- "RNA"
p15<-Downsample[[]] %>%
  subset(compartments == "Hematoma") %>%
  subset(Age_raw=="Immunoaged") %>%
  filter(CellTypes %in% c("B Cells")) %>% 
  select(sample,group, Age_raw,Genes_Bcells1,compartments,timepoint,CellTypes) %>% 
   #mutate(Age_raw=factor(Age_raw, levels = c("Young","Aged","ImmunoAged"))) %>%
  ggplot(aes(timepoint, Genes_Bcells1))+ 
  geom_boxplot(size = .4, outlier.size = 0) + 
  geom_point(aes(color = timepoint), size = .6) +
   #facet_wrap(~Age_raw, scales="free_x")+   
  ylim(0, 2) + 
  stat_compare_means(comparisons = my_comparisons,label="p.signif",
                  method = "t.test",label.y = 1.5, label.x = "5days",size = 6.1167,hide.ns = TRUE) + theme_classic() + 
theme(axis.title = element_text(size = 8), 
        axis.text = element_text(size = 8, family = "Arial"),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(),
        axis.text.y=element_text(size =8),
        legend.text = element_text(size =8),
        legend.title = element_text(size =8),
        strip.text.x = element_text(size = 8),
        strip.text.y = element_text(size = 8),
        plot.title = element_text(size =8,hjust = 0.5,))+
     scale_color_manual(name= "Time",values = c("dodgerblue3","orange2")) + labs(title ="Immune-aged: hematoma", x = c(color =''), y = 'Gene module score of LI.M47.0') 
p15

Generating Figure 3

(wrap_elements(p12)  + wrap_elements(p10)) /
  (wrap_elements(p11) + wrap_elements(p15) ) /
  (wrap_elements(p13) + wrap_elements(p14) ) +
  plot_layout(heights=c(1,.8,1.1)) +
  plot_annotation(tag_levels = 'A') &
  theme(plot.tag = element_text(size = 8, face = "bold",hjust = -.5, vjust = 0))

if (file.exists(file.path(figure_out_path,"fig3.pdf"))){
  try(ggsave(file.path(figure_out_path,"fig3.pdf"), device = cairo_pdf,width=8,height=9), silent=T)
  warning("file not written, please close file in the other application and rerun.")
} else {
  ggsave(file.path(figure_out_path,"fig3.pdf"), device = cairo_pdf,width=8,height=9)
}

Generating Suppl Figure 3

(wrap_elements(sf3a)  + wrap_elements(sf3b)) +
  plot_annotation(tag_levels = 'A') &
  theme(plot.tag = element_text(size = 8, face = "bold",hjust = -.5, vjust = 0))

if (file.exists(file.path(figure_out_path,"SFig3.pdf"))){
  try(ggsave(file.path(figure_out_path,"SFig3.pdf"), device = cairo_pdf,width=8,height=9), silent=T)
  warning("file not written, please close file in the other application and rerun.")
} else {
  ggsave(file.path(figure_out_path,"SFig3.pdf"), device = cairo_pdf,width=8,height=3)
}

supplement figure (SFig A to D)

# library(RColorBrewer)
# #################### Volcano plot #######################
# DF <- res_df[grepl("12w_nonSPF_5days_fx_vs_12w_nonSPF_2days_fx", res_df$contrast), ]
# (SFig4_1 <- DF %>% 
#     dplyr::filter(!is.na(padj)) %>%
#     ggplot(aes(x=log2FoldChange,y=-log10(pvalue),color=padj < .1)) + 
#     geom_point(size=.5, shape=1) +
#     geom_label_repel(data=DF %>%
#                        dplyr::filter(padj < .1),
#                      aes(label=gene_name),
#                      segment.size=.25,segment.alpha=.5,size=4,color='black',max.overlaps=15,
#                      label.padding=0.05,label.size=NA) +
#     scale_color_manual(values=c('black','red'))+ theme_classic() +  theme(text = element_text(family = "Arial"),
#         plot.title = element_text(size = 8,hjust = 0.5,),
#                      axis.title = element_text(size = 8),  
#                      axis.text = element_text(size = 8),
#                     axis.ticks.x=element_blank(),
#                      legend.text = element_text(size =8),
#                      legend.title = element_text(size = 8),
#                     legend.position = "none",
#                      strip.text.x = element_blank()) + labs(title = "Young",x = "log2FoldChange (5 days vs. 2 days)"))
# 
# DF <- res_df[grepl("26w_SPF_5days_fx_vs_26w_SPF_2days_fx", res_df$contrast), ]
# (SFig4_2 <- DF %>% 
#     dplyr::filter(!is.na(padj)) %>%
#     ggplot(aes(x=log2FoldChange,y=-log10(pvalue),color=padj < .1)) + 
#     geom_point(size=.5, shape=1) +
#     geom_label_repel(data=DF %>%
#                        dplyr::filter(padj < .1),
#                      aes(label=gene_name),
#                      segment.size=.25,segment.alpha=.5,size=4,color='black',max.overlaps=15,
#                      label.padding=0.05,label.size=NA) +
#     scale_color_manual(values=c('black','red'))+ theme_classic() +  theme(text = element_text(family = "Arial"),
#         plot.title = element_text(size = 8,hjust = 0.5,),
#                      axis.title = element_text(size = 8),  
#                      axis.text = element_text(size = 8),
#                     axis.ticks.x=element_blank(),
#                      legend.text = element_text(size =8),
#                      legend.title = element_text(size = 8),
#                     legend.position = "none",
#                      strip.text.x = element_blank()) + labs(title = "Aged",x = "log2FoldChange (5 days vs. 2 days)"))

Generating supplement Figure E

# library(tmod)
# library(disco)
# dataset_1 <- res_df[res_df$contrast == "12w_nonSPF_5days_fx_vs_12w_nonSPF_2days_fx", ]
# dataset_2 <- res_df[res_df$contrast == "26w_nonSPF_5days_fx_vs_26w_nonSPF_2days_fx", ]
# dataset_3 <- res_df[res_df$contrast == "26w_SPF_5days_fx_vs_26w_SPF_2days_fx", ]
# names<- c("Aged","Immunoaged")
# orthologs <- NULL
# orthologs <- makeMatchedOrtholog(names,dataset_3$gene_name, dataset_3$log2FoldChange,dataset_2$log2FoldChange, dataset_3$padj,dataset_2$padj,row.names = NULL, extra = NULL)
# #ds <- disco.score(orthologs)
# #plotDisco(orthologs, ds)
# 
# both<- orthologs$lfc %>% as.data.frame()
# both$gene_name <- orthologs$genes
# genes <- res_df[res_df$contrast == c("26w_SPF_5days_fx_vs_26w_SPF_2days_fx","26w_nonSPF_5days_fx_vs_26w_nonSPF_2days_fx"), ]
# Genes <- genes[genes$padj < 0.05  & 
#                               !is.na(genes$padj) &
#                               abs(genes$log2FoldChange) > 0.5,] %>% pull(gene_name)

Plot logfoldchange against one another

# # add a column of NAs
# both$colours <- "NO"
# # if log2Foldchange > 0.6 and pvalue < 0.05, set as "UP" 
# both$colours[(both[,1] > 0.5 & dataset_3$padj< 0.05)] <- "UP"
# both$colours[(both[,2] > 0.5 & dataset_2$padj< 0.05)] <- "UP"
# # if log2Foldchange < -0.6 and pvalue < 0.05, set as "DOWN"
# both$colours[(both[,1] < -0.5 & dataset_3$padj< 0.05)] <- "DOWN"
# both$colours[(both[,2] < -0.5 & dataset_2$padj< 0.05)] <- "DOWN"
# both$gene_set_labels <- ifelse(both$gene_name %in% Genes,both$gene_name, NA)
# SFig4_5 <- ggplot(both, aes(both[,1], both[,2],colour = colours)) + geom_point(shape=1, size = .5) + geom_text_repel(aes(label = gene_set_labels), size = 3, max.overlaps = 100) + ggtitle("Hematoma: Aged vs Immunoaged") + xlab("log2Foldchange (Aged)") + ylab("log2Foldchange (Immunoaged)") +  scale_color_manual(values=c("blue","grey","red")) + theme_classic() +  theme(plot.title = element_text(size = 8, hjust=0.5))
# SFig4_5

Supplemenyt Figure F: disco plot young vs immunoaged

# library(tmod)
# library(disco)
# dataset_1 <- res_df[res_df$contrast == "12w_nonSPF_5days_fx_vs_12w_nonSPF_2days_fx", ]
# dataset_2 <- res_df[res_df$contrast == "26w_nonSPF_5days_fx_vs_26w_nonSPF_2days_fx", ]
# dataset_3 <- res_df[res_df$contrast == "26w_SPF_5days_fx_vs_26w_SPF_2days_fx", ]
# names<- c("Young","Immunoaged")
# orthologs <- NULL
# orthologs <- makeMatchedOrtholog(names,dataset_1$gene_name, dataset_1$log2FoldChange,dataset_2$log2FoldChange, dataset_1$padj,dataset_2$padj,row.names = NULL, extra = NULL)
# #ds <- disco.score(orthologs)
# #plotDisco(orthologs, ds)
# 
# both<- orthologs$lfc %>% as.data.frame()
# both$gene_name <- orthologs$genes
# genes <- res_df[res_df$contrast == c("12w_nonSPF_5days_fx_vs_12w_nonSPF_2days_fx","26w_nonSPF_5days_fx_vs_26w_nonSPF_2days_fx"), ]
# Genes <- genes[genes$padj < 0.05  & 
#                               !is.na(genes$padj) &
#                               abs(genes$log2FoldChange) > 0.5,] %>% pull(gene_name)

Plot logfoldchange against one another

# # add a column of NAs
# both$colours <- "NO"
# # if log2Foldchange > 0.6 and pvalue < 0.05, set as "UP" 
# both$colours[(both[,1] > 0.5 & dataset_1$padj< 0.05)] <- "UP"
# both$colours[(both[,2] > 0.5 & dataset_2$padj< 0.05)] <- "UP"
# # if log2Foldchange < -0.6 and pvalue < 0.05, set as "DOWN"
# both$colours[(both[,1] < -0.5 & dataset_1$padj< 0.05)] <- "DOWN"
# both$colours[(both[,2] < -0.5 & dataset_2$padj< 0.05)] <- "DOWN"
# both$gene_set_labels <- ifelse(both$gene_name %in% Genes,both$gene_name, NA)
# SFig4_6 <- ggplot(both, aes(both[,1], both[,2],colour = colours)) + geom_point(shape=1, size = .5) + geom_text_repel(aes(label = gene_set_labels), size = 3, max.overlaps = 100) + ggtitle("Hematoma:Young vs Immunoaged") + xlab("log2Foldchange (Young)") + ylab("log2Foldchange (Immunoaged)") +  scale_color_manual(values=c("blue","grey","red")) + theme_classic() +  theme(plot.title = element_text(size = 8, hjust=0.5))
# SFig4_6

generating tmod plot for figure 4

# x <- genes$"26w_SPF_5days_fx_vs_26w_SPF_2days_fx"
# sgenes <- tmodDecideTests(g = x$gene_name, lfc = x$log2FoldChange, pval = x$padj, lfc.thr = 0.5, pval.thr = 0.1, mset = tmod_mm)
# 
# names(sgenes) <- "26w_SPF_5days_fx_vs_26w_SPF_2days_fx"
# p10 <- ggPanelplot(res.tmod.filtered, sgenes = sgenes,q_cutoff = 1e-12) + 
#   theme(text = element_text(family = "Arial"),
#         plot.title = element_text(size = 8,hjust = 0.5,),
#         axis.title = element_text(size = 8),
#         axis.text = element_text(size = 8),
#         axis.ticks.x=element_blank(),
#         legend.text = element_text(size =8),
#         legend.title = element_text(size = 8),
#         strip.text.x = element_blank()) + 
#   labs(title = "Aged: \nhematoma, \n5 days vs 2 days")
# p10

Generating supplement Figure 4

# library(patchwork)
# SFig4_1+ SFig4_2  + SFig4_6 + SFig4_5  +
#   plot_annotation(tag_levels = 'A') + plot_layout(ncol = 2)
# ggsave(file.path(figure_out_path,"SFig4.pdf"), device = cairo_pdf, width=8,height=9)
sessionInfo()
## R version 4.3.3 (2024-02-29)
## Platform: x86_64-conda-linux-gnu (64-bit)
## Running under: Ubuntu 22.04.3 LTS
## 
## Matrix products: default
## BLAS/LAPACK: /home/milekm/work/miniconda3/envs/monocle/lib/libopenblasp-r0.3.27.so;  LAPACK version 3.12.0
## 
## locale:
##  [1] LC_CTYPE=C.UTF-8       LC_NUMERIC=C           LC_TIME=C.UTF-8       
##  [4] LC_COLLATE=C.UTF-8     LC_MONETARY=C.UTF-8    LC_MESSAGES=C.UTF-8   
##  [7] LC_PAPER=C.UTF-8       LC_NAME=C              LC_ADDRESS=C          
## [10] LC_TELEPHONE=C         LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C   
## 
## time zone: Europe/Berlin
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] stats4    stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] ggpubr_0.6.0                tmod_0.50.13               
##  [3] patchwork_1.2.0             RColorBrewer_1.1-3         
##  [5] stringr_1.5.1               ggrepel_0.9.5              
##  [7] DESeq2_1.42.0               SummarizedExperiment_1.32.0
##  [9] Biobase_2.62.0              MatrixGenerics_1.14.0      
## [11] matrixStats_1.3.0           GenomicRanges_1.54.1       
## [13] GenomeInfoDb_1.38.1         IRanges_2.36.0             
## [15] S4Vectors_0.40.2            BiocGenerics_0.48.1        
## [17] cowplot_1.1.3               ggplot2_3.5.1              
## [19] gtools_3.9.5                tidyr_1.3.1                
## [21] dplyr_1.1.4                 Seurat_5.1.0               
## [23] SeuratObject_5.0.2          sp_2.1-4                   
## 
## loaded via a namespace (and not attached):
##   [1] RcppAnnoy_0.0.22        splines_4.3.3           later_1.3.2            
##   [4] plotwidgets_0.5.1       bitops_1.0-7            tibble_3.2.1           
##   [7] polyclip_1.10-7         XML_3.99-0.17           fastDummies_1.7.3      
##  [10] lifecycle_1.0.4         rstatix_0.7.2           globals_0.16.3         
##  [13] lattice_0.22-6          MASS_7.3-60.0.1         crosstalk_1.2.1        
##  [16] backports_1.5.0         magrittr_2.0.3          plotly_4.10.4          
##  [19] sass_0.4.9              rmarkdown_2.27          jquerylib_0.1.4        
##  [22] yaml_2.3.9              httpuv_1.6.15           sctransform_0.4.1      
##  [25] spam_2.10-0             spatstat.sparse_3.1-0   reticulate_1.38.0      
##  [28] pbapply_1.7-2           abind_1.4-5             zlibbioc_1.48.0        
##  [31] Rtsne_0.17              purrr_1.0.2             msigdbr_7.5.1          
##  [34] RCurl_1.98-1.16         GenomeInfoDbData_1.2.11 irlba_2.3.5.1          
##  [37] listenv_0.9.1           spatstat.utils_3.0-5    pheatmap_1.0.12        
##  [40] goftest_1.2-3           RSpectra_0.16-2         spatstat.random_3.3-1  
##  [43] fitdistrplus_1.2-1      parallelly_1.37.1       leiden_0.4.3.1         
##  [46] codetools_0.2-20        DelayedArray_0.28.0     DT_0.33                
##  [49] tidyselect_1.2.1        farver_2.1.2            spatstat.explore_3.3-1 
##  [52] jsonlite_1.8.8          progressr_0.14.0        ggridges_0.5.6         
##  [55] survival_3.7-0          tools_4.3.3             tagcloud_0.6           
##  [58] ica_1.0-3               Rcpp_1.0.13             glue_1.7.0             
##  [61] gridExtra_2.3           SparseArray_1.2.2       xfun_0.46              
##  [64] withr_3.0.0             fastmap_1.2.0           fansi_1.0.6            
##  [67] caTools_1.18.2          digest_0.6.36           R6_2.5.1               
##  [70] mime_0.12               colorspace_2.1-0        scattermore_1.2        
##  [73] tensor_1.5              spatstat.data_3.1-2     utf8_1.2.4             
##  [76] generics_0.1.3          data.table_1.15.4       httr_1.4.7             
##  [79] htmlwidgets_1.6.4       S4Arrays_1.2.0          uwot_0.2.2             
##  [82] pkgconfig_2.0.3         gtable_0.3.5            lmtest_0.9-40          
##  [85] XVector_0.42.0          htmltools_0.5.8.1       carData_3.0-5          
##  [88] dotCall64_1.1-1         scales_1.3.0            png_0.1-8              
##  [91] spatstat.univar_3.0-0   knitr_1.48              reshape2_1.4.4         
##  [94] nlme_3.1-165            cachem_1.1.0            zoo_1.8-12             
##  [97] KernSmooth_2.23-24      parallel_4.3.3          miniUI_0.1.1.1         
## [100] pillar_1.9.0            grid_4.3.3              vctrs_0.6.5            
## [103] RANN_2.6.1              gplots_3.1.3.1          promises_1.3.0         
## [106] car_3.1-2               xtable_1.8-4            cluster_2.1.6          
## [109] beeswarm_0.4.0          evaluate_0.24.0         cli_3.6.3              
## [112] locfit_1.5-9.9          compiler_4.3.3          rlang_1.1.4            
## [115] crayon_1.5.3            ggsignif_0.6.4          future.apply_1.11.2    
## [118] labeling_0.4.3          plyr_1.8.9              stringi_1.8.4          
## [121] viridisLite_0.4.2       deldir_2.0-4            BiocParallel_1.36.0    
## [124] babelgene_22.9          munsell_0.5.1           lazyeval_0.2.2         
## [127] spatstat.geom_3.3-2     Matrix_1.6-5            RcppHNSW_0.6.0         
## [130] future_1.33.2           shiny_1.8.1.1           highr_0.11             
## [133] ROCR_1.0-11             broom_1.0.6             igraph_2.0.3           
## [136] bslib_0.7.0